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Recently, Jouve et al [2] published the paper that presents the numerical benchmark for the 
solar dynamo models. Here, I would like to show a way how to get it with help of computer 
algebra system Maxima. This way was used in |4] to test some new ideas in the large-scale 
stellar dynamos. What you need are the latest version of Maxima-5.16.3 (preferable compiled 
against the fastest lisps like sbcl or cmucl-sse2) and some knowledge of the global (spectral) 
methods to solve the PDE eigenvalue problem. For the quite comprehensive introduction to 
these methods please look at the book by John Boyd [T]. The basic steps to solve the problem 
are: 

1. the mathematical formulation (equation+boundary conditions) 

2. choice the basis function and project equations to the basis 

3. find matrices (apply some integration procedure in case of Galerkin method) 

4. apply linear algebra 

The whole consideration is divided for two cases. As the first case we explore the largest 
free decay modes in the sphere which is submerged in vacuum. In this problem the all dynamo 
effects are neglected. As the second case I test the aQ dynamo in the solar convection zone 
with the tachocline included. 

Lets consider the spherical geometry. The evolution of the axisymmetric large-scale mag- 
netic field (LSMF), (B) = e^B + curl (where r is radius, 9- co-latitude, e^- the unit 

azimuthal vector) in the turbulent media subjected to the differential rotation in the spherical 
shell can be described with equations: 

_ i d{n,A) i f djrs^ _ d£_\ 

~dt ~ r d{r,e) ^ r ylh W J ' 

^ = rsinOSt (2) 
dt ' ^ ^ 

In equations above, the turbulent contribution is expressed through the components of the 
mean electromotive force (MEMF) £^ = (u' x b'), where u',b' are the small-scale fluctuated 
velocity and magnetic field respectively, Q = Q {r,6) - the given angular velocity distribution. 
For the sake of simplicity we restrict consideration to the case of aQ dynamo with isotropic 
turbulent diffusion. Hence, we have 

_ Vt d sin OB 
~ rsinO 09 

Se = (4) 

r or 

■ nc [9^^ sine d 1 dA\ ^„ . 



*Thc text were processed with Emaxima ( jhttp://Maxima.sf.net| 
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where rjT- turbulent diffusion, a (r, 6) - dimensionless function to model tlie a effect, G = 
Vlogp - stratification parameter. We adopt the model parameters given in [2]. The boundary 
drB 

conditions are — — = 0, A = - at the bottom and vacuum conditions - at the top of convec- 
or 

tion zone. For the computation all the equations are written in dimensionless form with new 
radial coordinate x = t/Rq. Moreover, to project equations on the basis function we use the 
coordinate transformation to the interval where basis is orthogonal. 

As the first step we consider the solutions for the largest free-decay modes. This intends 
to test the accuracy of the boundary conditions implementation procedure and the speed of 
convergence. We neglect all the dynamo effects in ([Tl [2| and return to equations: 



dB _ ld^{xB) 1 d 1 d {sin 9B) 

dt X dx"^ x"^ 36 sin 9 36 ' 
3 A _ 3'^ A sine 3 1 3 A 

~dt ~ Ih^^^mMlW' 

where for the sake of simplicity we have assumed that magnetic diffusivity is constant over the 
depths and equations are written in dimensionless form. Lets consider the integration domain 
on the radial coordinate to be x G [0, 1] and for latitude - from pole to pole. Having the 
regular conditions at the origin and vacuum boundary conditions at the top, we can get the 
analytical solutions of ([6| [?]). They can be expressed via the spherical Bessel functions. For the 
sake of simplicity we restrict ourselves to solutions for the largest modes. We decompose the 
eigenmodes as follows 

B{x,e,t) = e^*5^5^6-5f)(x)P^(cos(^)), (8) 

n m 

A{x,e,t) = e^*^^a"'"sin(^)^(l)(x)P^(cos(^)), (9) 

n m 

where are the associated Legendre polynomials and other functions were found via basis 
recombination of the Legendre polynomials, namely, 

= x(P2„+i(x)-Pi(x)), (10) 

(.) ^ . (.) - p. (.) "^"^^'p;:f,;^""^^'> ) . (11) 

Here, each element of basis satisfies the boundary conditions individually . For the largest 
decay modes we have S[^^ (x) cc ji (ax), A = — ~ 4.4934) and Sl"^^ (x) cc ji (irx), A = 
— vr^. Similar to [3j, the eigenvectors are scaled as S[^^ (.5) = 1 and S[^^ (1) = 1. The errors 
are measured as E (X) = \Xtrue — Xnum\ and E (B) = \'Btrue — ^numf dx. The results of 
benchmark are given in the Table [T| 



Next benchmark is related to one given in the paper [2]. The detail of the model can be 
found in that paper. We consider the results for test B, which is for aQ dynamo with external 
vacuum boundary conditions and jump othe magnetic diffusivity below the botom of convection 
zone. For the magnetic fields we have used the following set of the modes: 

\ ISn'^ + 2bn + 6 Ura^ + 2bn + 6 / 

oiA)f. _ ^ /^p I + 3) Pn+i jx) in' + 4n + 2m// + 1) P„+2 (x) \ 

""^ ^ V ^' + 4^ + 2m// + 4 n2 + 4n + 2m// + 4 J' ^ ' 
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E{B) 


E (A) ,A 


E{A) 


3 


3.83e-5 


6.849e-4 


1.08e-7 


3.98e-9 


4 


9.984e-8 


4.365e-9 


5.651e-ll 


1.255e-12 


5 


1.207e-10 


2.497e-12 


9.68e-14 


2.142e-16 


6 


8.707e-14 


9.068e-16 


7.72e-14 


2.136e-20 


7 


1.77e-14 


4.80e-19 


1.38e-14 


1.325e-24 


8 


5.32e-15 


6.263e-23 


1.90e-14 


4.427e-29 



Table 1: Convergence of the decay modes to the larges mode of analytic solution of 
the number of modes in radial basis. 



Resolution 


^crit 


UJ 


8x8 


.443 


180.5 


10x10 


.4175 


175.1 


12x12 


.4095 


172.2 


13x13 


.411 


172.4 


14x14 


.4122 


172.7 


16x16 


.4125 


172.9 



Table 2: Model B of [2]. Resolution N x M means that N modes is used in decomposition on 
radius and M - on latitude. 

2 

where factor / = appears due the coordinate transformation [xi,l] [—1,1] and 

1 

Xi = 0.65 is bottom of integration domain. 

The results are shown at the table 2 and Figure [1} Note, that our results are quite close to 
those in [2|- The small difference can be attributed to the difference in the method of solution 
of the problem. 

The largest decay mode of the diffusion operator 

Here, I give some details of the above computation within Maxima. The strongest side of 
the computer algebra system is that most of the computational work is done exactly. For 
example, all the derivatives in (|6| [T]) are computed directly by substitution of ([8| [9]) to ([6} 
[T]). The integration over space domain has to be done to proceed with the Galerkin method. 
This is carried out with help of the Gauss-Legendre procedure. The procedure is based on 
the so-called "collocation points and weights" method. The collocation points are the set of 
zeros of the Legendre polynomial of the order > where is the number of the greatest 
mode among n = 1, . . . N in (jsj |9]). The collocation points and weights can be found via 
subprogram "pseudp.mac" . It is given in Appendix. The maxima session is started with the 
general definitions: 

kill (all) $ 

batchload ( "pseudp . mac " ) $ 

(xb:0. ,xe: 1.)$ 

bf torat :true$ 

f loat2bf :true$ 

ratprint : f alse$ 

fpprec:64$ 

ratepsilon : 1 . e-32$ 



3 



Figure 1: Model B of [2]. Snapshots of evolution of torioidal magnetic fileds (shadowed density 
plot) and poloidal fields (shown with streamlines) over the half of cycle. 

chn(x) : =x$ 

realGr(Ll ,L2) :=(if realpart(Ll) >= realpart(L2) then true else false)$ 
Here, we find the collocation points and weights on radius, 

(np : 5 , Nch : 4*np , xchwh : gaulegP (-1,1, Nch , 1 ) , xch : xchwh [1] , wh : xchwh [2] ) $ 
and on latitude, 

(NT : 6 , NTT : 6*NT , mut.w : gaulegP (-1,1, NTT , 1 ) , mut : mut_w [1] , w : mut.w [2] ) $ 

Next, we define the indices for the basis functions and the parities for the toroidal and poloidal 
eigenmodes, 

Nch: length (xch) $ 
NTT:length(mut)$ 

for mil thru NT do( for n:l thru np do( 

ki:np*(ni-l)+n,ni [ki] :n,mi [ki] :m))$ 
ni : makelist (ni [i] , i , 1 , NT*np) $ 
mi : makelist (mi [i] , i , 1 , NT*np) $ 
NN : NT*np$ 
parity_b:0$ 
parity_a: 0$ 
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Note, that only a half of [—1, 1] intervals is used in computations. This leads to the computa- 
tional economy and to increase the spectral resolution of the code. The following are the part 
of the code which defines the basis functions on radius {Sn^^ — > chtbn, 5'nm —>■ chtanm ) and 
latitude, 

kill(chtbl,chtb)$ 

chtb [n] (x) : =x* (legendre.p ( 1 , x) -legendre_p (2*n+l , x) ) $ 
makelist (taylor (chtb [i] (x) ,x,0,l) ,i,l,np)$ 
plot2d(makelist(chtb[i] (x) ,i,l,np) , [x,-l,l])$ 
kill(chtaOO)$ 

chtaOO(n,l,x) : =x* (legendre_p(2*n+l ,x) 

- ( (2*n+l) * (2*n+2) +2*1+2) *legendre_p (1 , x) / (2*1+4) ) $ 
chta[n,l] (x) :=chtaOO(n, 2*1-1, x) , 
kill(P110,Pla,Plb)$ 

PllO (n, x) : =expand( (sqrt ( (2*n+l) / (2*n* (n+1) ) ) *assoc_legendre_p(n, 1 ,x) ) ) $ 

Pla[n] (x) :=block(if parity_a=0 then bf loat (P110(2*n-1 ,x) ) else bf loat(P110(2*n,x)) ) $ 

Plb[n] (x) :=block(if parity_b=0 then bfloat(P110(2*n-l,x)) else bf loat(P110(2*n,x)) ) $ 

As you see, we use odd or even associated Legendre polynomial in respective of the parity 
choice. To accelerate the integration procedure we calculate the matrices of basis functions 
over the sets of collocation points, 

remarray (PLB , PLA) $ 
kill(PLBl,PLAl)$ 

PLB[n,i] :=bfloat(Plb[n] (mut[i]))$ 
PLA[n,i] :=bfloat(Pla[n] (mut[i]))$ 
PLAl [n , i] : =bf loat (Pla [n] (mut [i] ) *w [i] ) $ 
PLBl [n , i] : =bf loat (Plb [n] (mut [i] ) *w [i] ) $ 
CHTB[n, j] : =bf loat (chtb [n] (xch[j] ))$ 
CHTA[n,ni,j] :=bf loat (chta[n,m] (xch[j]))$ 

makelist(makelist(makelist(CHTA[n,m, j] ,n,l,np) ,m,l,NT) , j ,l,Nch)$ 
CHTBl[n,j] : =bf loat (chtb [n] (xch [j] ) *wh [j] ) $ 
CHTAl [n,m, j] :=bf loat(chta[n,m] (xch[j] )*wh[j] )$ 

makelist (makelist (makelist (CHTAl [n,in,j] ,n,l,np) ,m,l,NT) , j ,l,Nch)$ 

genmatr ix (PLB , NT , NTT) $ 

genmatr ix (PLA , NT , NTT) $ 

genmatr ix (PLAl , NT , NTT) $ 

genmatr ix (PLBl , NT , NTT) $ 

genmatr ix (PLB2 , NT , NTT) $ 

genmatr ix (CHTB , np , Nch) $ 

genmatr ix (CHTBl , np , Nch) $ 

Now we can to proceed to solution of the problem ([6} [?]). They can be treated separately. 
Firstly we define the left parts of 



kill(AA,BB)$ 

AA[ki,kj] :=bfloat(sum(CHTAl[ni [ki] ,mi[ki] ,j]*CHTA[ni [kj] ,mi[kj] ,j]* 

sum (PLAl [mi [ki] , i] *PLA [mi [kj] ,i] ,i,l,NTT) ,j ,l,Nch))$ 
BB[ki,kj] :=bfloat(sum(CHTBl[ni [ki] , j] *CHTB [ni [kj] ,j] 
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*sum(PLBl[mi [ki] ,i]*PLB[mi [kj] ,i] , i , 1 ,NTT) , j , 1 ,Nch) )$ 
BB : genmatrix (BB , NN , NN) $ 
AA : genmatrix (AA , NN , NN) $ 

Note, the integration is just a product of sums over collocation points. Now lets compute the 
right part of ([6|. The computations are straightforward. The radial part is, 



kill(CHTB_d2)$ 

CHTB_d2[n,j] : =bf loat (subst ( [x=xch [j] ] , 

expand((diff (diff (chtb[n] (x) ,x) ,x) ))))$ 
genmatrix (CHTB_d2 , np , Nch) $ 
kill(et_d)$ 

et_d[ki,kj] :=bf loat (sum (CHTB_d2 [ni [kj] , j] *CHTB1 [ni [ki] , j] , j ,l,Nch) 

*sum(PLB[mi [kj] ,i]*PLBl[mi [ki] ,i] ,i,l,NTT))$ 
et_d : genmatrix (et_d , NN , NN) $ 

The latitudinal part is computed as follows 



kill (PLB_d2 , PLB_d2m , PLB.dl , PLB.d) $ 

PLB_d2m[n,i] : =bf loat (subst ( [mu=mut [i] ] , expand( (sqrt (l-mu~2) 

*diff (diff (sqrt(l-mu"2)*Plb [n] (mu) ,mu) ,mu) ) ) ) ) $ 
genmatrix (PLB_d2m , NT , NTT) $ 
kill(er_d)$ 

er_d [ki , kj ] : =bf loat (sum (CHTBl [ni [ki] , j ] *CHTB [ni [kj ] , j ] /chn (xch [ j ] ) "2 , j , 1 , Nch) 

*sum(PLBl[mi [ki] , i] *PLB_d2m [mi [kj] ,i] ,i,l,NTT))$ 
er_d : genmatrix (er_d , NN , NN) $ 

Invert the left part of and solve the eigenvalue problem with help of lapack 



BMI : invert_by_lu (BB) $ 

realGrO(Ll) :=(if realpart(Ll) >=0 then true else false)$ 
load(lapack)$ 
MT:BMI. (er_d+et_d)$ 
lamb: (dgeev(MT,true) )$ 



; ; Loading #P"/usr/share/maxima/5 . 16 . 3/share/lapack/lapack/binary-cmucl/dtrevc . sse2f " . 

Now, lets sort the eigenvalues, 

lmb:sort(flatten((col(lamb,l)) [1]) ,realGr) ; 



[- 20.19072855751148, - 48.83125052872499, - 59.68036461428608, 

- 87.54250816953108, - 108.8085726200047, - 119.8731393320569, 

- 135.9006458324267, - 169.2937978826357, - 193.7765731289586, 

- 203.8596788530381, - 238.1199364033153, - 241.2729384729356, 

- 260.8709651925637, - 277.930972796181, - 320.6253098224506, 
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- 328.0254017082749, - 373.0041226108443, - 410.8018157235754, 

- 505.412713312208, - 577.4557504236772, - 665.4381736491396, 

- 809.3382250971157, - 851.3024977594987, - 882.9750962899574, 

- 932.6139588501989, - 970.2769707279747, - 1049.716092377411, 

- 1259.650734653339, - 1830.07144488619, - 2585.672953057317] 

The relevant spherical Bessel functions of the problem are 

jn(n,x) : =spherical_bessel_j (n,x)$ 
jnR(n,x) : =x*spherical_bessel_j (n,x)$ 

Next we compute the first root of jn{x) and compare it with the first eigenvalue, 

al:f ind_root(jn(l,x) ,x,4,5)$ 
sqrt (abs (f irst (Imb) ) ) -al ; 



1 . 207158817351228e-10 
Similarly, we solve ([T]). At the first, we invert the left part of it 

AMI : invert_by_lu (AA) $ 
Then calculate the right part and apply lapack solver, 

kill (CHTA_d2 , PLA_d2) $ 

CHTA_d2[n,in, j] : =subst ( [x=xch [j] ] , (diff (diff (chta[n,m] (x) ,x) ,x)))$ 
makelist(makelist(makelist(CHTA_d2[n,ni, j] ,n,l,np) ,m,l,NT) , j ,l,Nch)$ 
PLA_d2[n,i] :=bf loat(subst( [mu=mut [i]] , (sqrt (l-mu"2) 
*diff (diff (sqrt (l-mu"2)*Pla[n] (mu) ,mu) ,mu))))$ 
genmatr ix (PLA_d2 , NT , NTT) $ 
kill(ef_d)$ 

ef_d[ki,kj] :=bf loat (sum(CHTA_d2 [ni [kj] ,mi[kj] ,j] 
*CHTA1 [ni [ki] , mi [ki] , j ] , j , 1 , Nch) * 
suin(PLAl[mi [ki] ,i]*PLA[nii [kj] ,i] ,i,l,NTT) 
+sum(CHTA[ni [kj] ,mi[kj] , j] *CHTA1 [ni [ki] ,nii[ki] ,j] 
/chn (xch [ j ] ) ~2 , j , 1 , Nch) * 

sum(PLAl[ini[ki] , i] *PLA_d2 [mi [kj] , i] ,i,l,NTT))$ 
ef _d : genmatrix (ef _d , NN , NN) $ 
MT: (AMI.ef_d)$ 
lamb: (dgeev(MT,true) )$ 

Then, we sort the eigenvalues and compare the largest eigen mode with the first zero of Ji/2{x), 

lmbl:sort(flatten((col(lamb,l)) [1] ) ,realGr)$ 
float (sqrt (abs (first (Imbl) ) ) -"/opi) ; 



- 9.681144774731365e-14 

To compare the eigenfunctions we have to apply the scaling s[^^ (.5) = 1 and S[^^ (1) = 1. 

The solution of dynamo problem is obtained in essentially the same way. I will put the code 
for it somewhere in the public place. 
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Appendix 

The subroutine to find the collocation points and weghts for the polynomial order of n in 
interval [xl,x2]. The Newton method is apphed, 

load(orthopoly)$ 

orthopoly_returns_intervals : f alse$ 

fpprec:32$ 

float2bf :true$ 

ratprint : f alse$ 

load ( " linearalgebra" ) $ 

gaulegP(xl,x2,n) := block( [xm,wm,AO,i,inm,pp,xin,xl,z,zl,numer,eps] , 

kill (xt , w) , eps : 1 . e-16 , f pprec : 32 , 

(if oddp(n) then mm:(n+l)/2 else min:n/2), 

xm:0.5*(x2+xl) , xl : . 5* (x2-xl) , 

for i: 1 while i <= mm do ( 

z : bf loat (cos (%pi* (i-0 . 25) / (n+0 . 5) ) ) , 

do (pp:bf loat (n*(z*legendre_p(n,z)-legendre_p(n-l ,z) )/(z*z-l . 0) ) , 

zl : z,z :bf loat (zl-legendre_p(n,z)/pp) , 

if abs(z-zl) < eps 

then return(zl)), 

xt [i] :xm-xl*zl, 

xt [n+l-i] :xm+xl*zl) , 

A :bf loat (genmatrix (lambda ( [i , j] ,legendre_p(i-l ,xt [j] ) ) ,n,n) ) , 
RH : bf loat (transpose (matrix ( 

makelist((if i=l then def int (legendre_p(0 ,x) ,x,-l , 1) else 0) , i , 1 ,n) ) ) ) , 
w: bf loat (invert _by_lu (A) . RH) , 

return ( [makelist (xt [i] , i , 1 ,n) , flatten (makelist(w[i] ,i,l,n))]))$ 
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